#ifndef ESMD_TYPES_H_
#define ESMD_TYPES_H_
/*defined type for use*/
typedef double real;
typedef long realint;
#include <limits.h>
typedef long tagint;
constexpr long TAG_MIN = LONG_MIN;
constexpr long TAG_MAX = LONG_MAX;
#include "dimension.h"
#define mpi_real MPI_DOUBLE
#ifndef __always_inline
#define __always_inline __attribute__((always_inline)) inline
#endif
bool is_finite_real(real x){
  long bits = *(long*)&x;
  return (bits >> 52 & 0x7ff) < 0x7ff;
}
/////////////////////////////////////////////////////////////////////////
// #ifndef __always_inline                                             //
// #define __always_inline __attribute__((artificial)) __always_inline //
// #endif                                                              //
/////////////////////////////////////////////////////////////////////////
#include <cmath>
#include <type_traits>
// #define ESMD_SP
// #ifdef ESMD_SP
// typedef float real;
// #define mpi_real MPI_FLOAT
// #ifdef __sw__
// #ifdef __sw_slave__
// #include <simd.h>
// typedef floatv4 vreal;
// typedef int realint;
// #define simd_vsqrt simd_vsqrts
// #endif
// #endif
// #else
// typedef double real;
// #define mpi_real MPI_DOUBLE
// #ifdef __sw__
// #ifdef __sw_slave__
// #ifdef __sw5__
// #include <simd.h>
// typedef doublev4 vreal;
// typedef long realint;
// #define simd_vsqrt simd_vsqrtd
// #else
// // #ifdef __sw7__
// // #include <simd.h>
// // typedef double vreal __attribute__ ((__mode__(__V8DF__)));
// // typedef long realint;
// // #define simd_vsqrt simd_vsqrtd
// // #endif
// #endif
// #endif
// #endif
// #endif
//typedef double double;
#define F_LR "%lf"
#define F_GR "%lf"
#define F_LI "%d"
#define F_GI "%ld"
#define REP2(x) x x
#define REP3(x) x x x
#define REP4(x) x x x x
#define REP5(x) x x x x x
#define REP6(x) x x x x x x

#define MAX_ELEM_NAME_LEN 7
typedef char elemstr[MAX_ELEM_NAME_LEN + 1];
//e1str-e2str-e3str-e4str'\0'
typedef char topstr[MAX_ELEM_NAME_LEN * 4 + 4];
typedef int typeint;
//typedef unsigned long unsigned long;
typedef long bond_t[2], angle_t[3], dihed_t[4];
#define CAST(type, input) ((type)(input))

/*vector type*/
template <typename T1, typename T2>
using OpType = decltype((T1)0 + (T2)(0));
template<typename T>
__always_inline T tsat(T x) {
  if (x < 0) return -x; else return x;
}
template <typename T>
struct vec{
  T x, y, z;
  vec(){};
  vec(T x, T y, T z) : x(x), y(y), z(z) {};
  template<typename OT> __always_inline vec<OpType<T, OT>> operator+(const vec<OT> &other) const {
    return vec<OpType<T, OT>>(x + other.x, y + other.y, z + other.z);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator-(const vec<OT> &other) const {
    return vec<OpType<T, OT>>(x - other.x, y - other.y, z - other.z);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator*(const vec<OT> &other) const {
    return vec<OpType<T, OT>>(x * other.x, y * other.y, z * other.z);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator/(const vec<OT> &other) const {
    return vec<OpType<T, OT>>(x / other.x, y / other.y, z / other.z);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator+(const OT &other) const {
    return vec<OpType<T, OT>>(x + other, y + other, z + other);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator-(const OT &other) const {
    return vec<OpType<T, OT>>(x - other, y - other, z - other);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator*(const OT &other) const {
    return vec<OpType<T, OT>>(x * other, y * other, z * other);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator/(const OT &other) const {
    return vec<OpType<T, OT>>(x / other, y / other, z / other);
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> cross(const vec<OT> &other) const {
    return vec<OpType<T, OT>>(y * other.z - other.y * z, other.x * z - x * other.z, x * other.y - other.x * y);
  }
  template<typename OT> __always_inline OpType<T, OT> dot(const vec<OT> &other) const {
    return x * other.x + y * other.y + z * other.z;
  }
  template<typename OT> __always_inline vec<T> &operator=(const vec<OT> &other) {
    x = other.x;
    y = other.y;
    z = other.z;
    return *this;
  }
  template<typename OT> __always_inline vec<T> &operator=(const OT &other) {
    x = other;
    y = other;
    z = other;
    return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator+=(const vec<OT> &other) {
    set(x + other.x, y + other.y, z + other.z); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator-=(const vec<OT> &other) {
    set(x - other.x, y - other.y, z - other.z); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator*=(const vec<OT> &other) {
    set(x * other.x, y * other.y, z * other.z); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator/=(const vec<OT> &other) {
    set(x / other.x, y / other.y, z / other.z); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator+=(const OT &other) {
    set(x + other, y + other, z + other); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator-=(const OT &other) {
    set(x - other, y - other, z - other); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator*=(const OT &other) {
    set(x * other, y * other, z * other); return *this;
  }
  template<typename OT> __always_inline vec<OpType<T, OT>> operator/=(const OT &other) {
    set(x / other, y / other, z / other); return *this;
  }
  __always_inline vec<T> operator-(){
    return vec<T>(-x, -y, -z);
  }
  template<typename OT> __always_inline void set(const OT &ox, const OT &oy, const OT &oz) {
    x = ox;
    y = oy;
    z = oz;
  }
  template<typename FT> __always_inline vec<FT> map(FT (*f)(T)) const {
    return vec<FT>(f(x), f(y), f(z));
  }
  __always_inline T vol() const {
    return x * y * z;
  }
  __always_inline T norm2() const {
    return x * x + y * y + z * z;
  }
  __always_inline T norm() const {
    return sqrt(norm2());
  }
  __always_inline vec<T> floor() const {
    return map<T>(std::floor);
  }
  __always_inline vec<T> unit_vec() const {
    return (*this) / this->norm();
  }
  vec<T> abs() const {
    return map<T>(std::abs);
  }
  vec<T> sat() const {
    return map<T>(tsat<T>);
  }
};
template<typename T, typename OT> __always_inline vec<OpType<T, OT>> operator+(const OT &s, const vec<T> &v) {
    return vec<OpType<T, OT>>(s + v.x, s + v.y, s + v.z);
}
template<typename T, typename OT> __always_inline vec<OpType<T, OT>> operator-(const OT &s, const vec<T> &v) {
    return vec<OpType<T, OT>>(s - v.x, s - v.y, s - v.z);
}
template<typename T, typename OT> __always_inline vec<OpType<T, OT>> operator*(const OT &s, const vec<T> &v) {
    return vec<OpType<T, OT>>(s * v.x, s * v.y, s * v.z);
}
template<typename T, typename OT> __always_inline vec<OpType<T, OT>> operator/(const OT &s, const vec<T> &v) {
    return vec<OpType<T, OT>>(s / v.x, s / v.y, s / v.z);
}
// #define defvec(base)  \
//   struct vec_##base { \
//     base x, y, z;     \
//   }

// #define vec(base) struct vec_##base
// #define vec(base) vec<base>
// #define box(base) Box<base>
// #define defbox(base)  \
//   struct box_##base { \
//     vec(base) lo, hi; \
//   }

// #define box(base) struct box_##base
#define vecarr(v0) \
  { (v0).x, (v0).y, (v0).z }
#define veccpy(v0, v1) (v0) = (v1)
#define vecfloor(v0, v1) (v0) = (v1).floor()
#define vecscale(v0, v1, scale) (v0) = (v1) * (scale)
#define vecadd(v0, v1, add) (v0) = (v1) + (add)
#define vecscaleadd(v0, v1, scale, add) (v0) = (v1) * (scale) + (add)

#define vecabs(v0, v1) (v0) = (v1).abs()
#define vecsat(v0, v1) (v0) = (v1).sat()
#define vecaddv(v0, v1, v2) (v0) = (v1) + (v2)
#define vecsubv(v0, v1, v2) (v0) = (v1) - (v2)
#define vecdivv(v0, v1, v2) (v0) = (v1) / (v2)
#define vecmulv(v0, v1, v2) (v0) = (v1) * (v2)
#define vecmuladdv(v0, v1, v2, v3) (v0) = (v1) * (v2) + (v3)
#define vecscaleaddv(v0, v1, v2, s1, s2) (v0) = (v1) * (s1) + (v2) * (s2)
#define vecset1(v, s) v.set(s, s, s)
#define vecset3(v, sx, sy, sz) v.set(sx, sy, sz)
#define vecvol(v) v.vol() 
#define vecnorm2(v) v.norm2() //((v).x * (v).x + (v).y * (v).y + (v).z * (v).z)
#define vecnorm(v) v.norm() //(sqrt(vecnorm2(v)))
#define vecdot(v0, v1) (v0).dot(v1) //((v0).x * (v1).x + (v0).y * (v1).y + (v0).z * (v1).z)
#define veccross(v0, v1, v2) (v0) = (v1).cross(v2)
#define veccross3(v0, v1, v2, v3) (v0) = (v1).cross(v2).cross(v3)

#define vecmaxv(v0, v1, v2)       \
  {                               \
    (v0).x = max((v1).x, (v2).x); \
    (v0).y = max((v1).y, (v2).y); \
    (v0).z = max((v1).z, (v2).z); \
  }
#define vecminv(v0, v1, v2)       \
  {                               \
    (v0).x = min((v1).x, (v2).x); \
    (v0).y = min((v1).y, (v2).y); \
    (v0).z = min((v1).z, (v2).z); \
  }
#define vecmax(v0, v1, scala)    \
  {                              \
    (v0).x = max((v1).x, scala); \
    (v0).y = max((v1).y, scala); \
    (v0).z = max((v1).z, scala); \
  }
#define vecmin(v0, v1, scala)    \
  {                              \
    (v0).x = min((v1).x, scala); \
    (v0).y = min((v1).y, scala); \
    (v0).z = min((v1).z, scala); \
  }

template<typename T>
struct box {
  vec<T> lo, hi;
  box(){}
  box(const vec<T> &lo, const vec<T> &hi): lo(lo), hi(hi) {}
  template<typename OT> __always_inline box<OpType<T, OT>> operator+(const OT &other) const {
    return box<OpType<T, OT>>(lo + other, hi + other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator-(const OT &other) const {
    return box<OpType<T, OT>>(lo - other, hi - other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator*(const OT &other) const {
    return box<OpType<T, OT>>(lo * other, hi * other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator/(const OT &other) const {
    return box<OpType<T, OT>>(lo / other, hi / other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator+(const vec<OT> &other) const {
    return box<OpType<T, OT>>(lo + other, hi + other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator-(const vec<OT> &other) const {
    return box<OpType<T, OT>>(lo - other, hi - other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator*(const vec<OT> &other) const {
    return box<OpType<T, OT>>(lo * other, hi * other);
  }
  template<typename OT> __always_inline box<OpType<T, OT>> operator/(const vec<OT> &other) const {
    return box<OpType<T, OT>>(lo / other, hi / other);
  }

  __always_inline T vol() const {
    return (hi - lo).vol();
  }
  template<typename OT>
  __always_inline bool contains(vec<OT> v) {
    return (v.x < hi.x && v.x >= lo.x && v.y < hi.y && v.y >= lo.y && v.z < hi.z && v.z >= lo.z);
  }
  template<typename OT>
  __always_inline box<T> &operator=(const box<OT> &other) const {
    lo = other.lo;
    hi = other.hi;
  }
};
#define boxcpy(b0, b1) (b0) = (b1)
#define boxscale(b0, b1, scale) (b0) = (b1) * (scale)
#define boxadd(b0, b1, add) (b0) = (b1) + (add)
#define boxscaleadd(b0, b1, scale, add) (b0) = (b1) * (scale) + (add)
#define vec_in_box(v, b) (b).contains(v)

// defvec<real>;
// defvec<int>;
// defvec<double>;
// defvec<float>;
template<typename T0, typename T1>
__always_inline OpType<T0, T1> min(T0 a, T1 b) {
  return ((a) < (b) ? (a) : (b));
}
template<typename T0, typename T1>
__always_inline OpType<T0, T1> max(T0 a, T1 b) {
  return ((a) > (b) ? (a) : (b));
}
// #define max(a, b) ((a) > (b) ? (a) : (b))
// #define min(a, b) ((a) < (b) ? (a) : (b))
#define max3(a, b, c) (max(max(a, b), c))
#define max4(a, b, c, d) (max(max(a, b), max(c, d)))
#define min3(a, b, c) (min(min(a, b), c))
#define min4(a, b, c, d) (min(min(a, b), min(c, d)))

typedef struct mdstat {
  double ecoul, evdwl, ebond, eangle, etori, eimpr;
  double virial[6];
  double need_rebuild;
} mdstat_t;
#define M_NA 6.02214076e23
#define M_K 1.38064852e-23
#include <math.h>
#ifndef M_PI
#define M_PI 3.141592653
#endif
#endif
